Introduction

In this script we are going to compute marker genes for the samples at varying resolutions.

Libraries

library(Seurat)
library(ggpubr)
library(cowplot)
library(dplyr)
library(ggplot2)
library(stringr)
library(readr)

Setting parameters

Loading necessary paths and parameters

source(here::here("misc/paths.R"))
source(here::here("utils/bin.R"))

set.seed(123)
source(here::here("misc/paths.R"))
source(here::here("utils/bin.R"))

"{anot}/{plt_dir}" %>%
  glue::glue() %>%
  here::here() %>%
  dir.create(path = .,
             showWarnings = FALSE,
             recursive = TRUE)

"{anot}/{robj_dir}" %>%
  glue::glue() %>%
  here::here() %>%
  dir.create(path = .,
             showWarnings = FALSE,
             recursive = TRUE)

Load data

The data used in this Rmarkdown document comes from 03-clustering_integration.Rmd where the data was integrated.

merged_se <- "{clust}/{robj_dir}/integrated_spatial.rds" %>%
  glue::glue() %>%
  here::here() %>%
  readRDS(file = .)

Analysis

First we are going to look at the cluster resolution we are going to check

Seurat::DimPlot(merged_se, group.by = "Spatial_snn_res.0.3")

res_vec <- c("Spatial_snn_res.0.01", "Spatial_snn_res.0.05", 
             "Spatial_snn_res.0.1", "Spatial_snn_res.0.3",
             "Spatial_snn_res.0.8", "Spatial_snn_res.1",
             "Spatial_snn_res.1.25", "Spatial_snn_res.1.5")
# Iterate over resolutions
lapply(res_vec, function(i) {
  # Iterate over slices
  lapply(id_sp_df$gem_id, function(id) {
    Seurat::SpatialPlot(
      merged_se,
      group.by = i,
      images = id,
      image.alpha = 0,
      crop = FALSE,
      pt.size.factor = 1.25)
  }) %>% 
    patchwork::wrap_plots(., nrow = 2, guides = "collect")
})
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

## 
## [[8]]

We are going to use the function FindAllMarkers to get the marker genes for all the clusters at the varying resolutions:

# Assign Identities
Seurat::Idents(merged_se) <- merged_se@meta.data[, "Spatial_snn_res.0.3"]

# Compute markers
markers03 <- Seurat::FindAllMarkers(
  object = merged_se,
  assay = "Spatial",
  slot = "data",
  only.pos = TRUE)

DT::datatable(
  data = markers03,
  filter = "top")

Save markers for resolution 0.3

out_file03 <- "{anot}/{robj_dir}/markers_Spatial_snn_res.0.3.xlsx" %>%
  glue::glue() %>%
  here::here()

# Remove file to write it from scratch
file.remove(out_file03)

lapply(unique(as.character(markers03$cluster)), function(clust) {
  
  markers03 %>%
    dplyr::filter(cluster == clust) %>%
    dplyr::arrange(dplyr::desc(avg_log2FC)) %>%
    xlsx::write.xlsx(
    .,
    file = out_file03,
    row.names = FALSE,
    sheetName = glue::glue("Cluster-{clust}"),
    append = file.exists(out_file03))
  })

Compute markers Resolution 1.5

We are going to use the function FindAllMarkers to get the marker genes for all the clusters at the varying resolutions:

# Assign Identities
Seurat::Idents(merged_se) <- merged_se@meta.data[, "Spatial_snn_res.1.5"]

# Compute markers
markers_15 <- Seurat::FindAllMarkers(
  object = merged_se,
  assay = "Spatial",
  slot = "data",
  only.pos = TRUE)

DT::datatable(
  data = markers_15,
  filter = "top")

Save markers for resolution 1.5

out_file_15 <- "{anot}/{robj_dir}/markers_Spatial_snn_res.1.5.xlsx" %>%
  glue::glue() %>%
  here::here()

# Remove file to write it from scratch
file.remove(out_file_15)

lapply(unique(as.character(markers_15$cluster)), function(clust) {
  markers_15 %>%
    dplyr::filter(cluster == clust) %>%
    dplyr::arrange(dplyr::desc(avg_log2FC)) %>%
    xlsx::write.xlsx(
    .,
    file = out_file_15,
    row.names = FALSE,
    sheetName = glue::glue("Cluster-{clust}"),
    append = file.exists(out_file_15))
  })

Clsutering Annotation

merged_se[["annotation-general"]] <- dplyr::case_when(
  merged_se@meta.data$Spatial_snn_res.0.3 == 0 ~ "T cell zone",
  merged_se@meta.data$Spatial_snn_res.0.3 == 1 ~ "T:B Border",
  merged_se@meta.data$Spatial_snn_res.0.3 == 2 ~ "Epithelial-1",
  merged_se@meta.data$Spatial_snn_res.0.3 == 3 ~ "Inter-Follicular Zone",
  merged_se@meta.data$Spatial_snn_res.0.3 == 4 ~ "Proliferating Follicle",
  merged_se@meta.data$Spatial_snn_res.0.3 == 5 ~ "Germinal Center",
  merged_se@meta.data$Spatial_snn_res.0.3 == 6 ~ "Subepithelial",
  merged_se@meta.data$Spatial_snn_res.0.3 == 7 ~ "Epithelial-2",
  merged_se@meta.data$Spatial_snn_res.0.3 == 8 ~ "Muscle"
)

Histological Annotation

Furthermore, we annotate these slides following a pathologists histological annotation. As of now we only have annotation of specific regions from one sample:

# load annotations
ann_dir <- "{anot}/{robj_dir}/mantle-zone-and-vascular-spaces" %>%
  glue::glue() %>%
  here::here()

hist_df <- lapply(list.files(ann_dir), function(n) {
  id <- substr(n, 1, 13)  # Extract first 13 characters
  
  # Go into the directory for each slide
  z_ls <- list.files(glue::glue("{ann_dir}/{n}"), full.names = TRUE)
  tmp_df <- lapply(z_ls, function(z) {
    
    # Iterate over region directories
    fn_ls <- list.files(z, full.names = TRUE)
    z_bc <- lapply(fn_ls, function(fn)
      readr::read_csv(file = fn) %>% dplyr::pull("barcode")) %>%
      unlist()
    
    # Add region to data frame
    dirs <- stringr::str_split(string = z, pattern = "/", simplify = TRUE)
    zone <- dirs[, ncol(dirs)]
    data.frame("barcodes" = z_bc, "histological_annotation" = zone)
  }) %>% dplyr::bind_rows()
  # Add gem_id to dataframe
  tmp_df$gem_id <- id
  tmp_df
  }) %>%
  dplyr::bind_rows()

# Make sure there are no duplicate rows
hist_df <- dplyr::distinct(hist_df)

Add histological annotation to seurat object

merged_se@meta.data <- merged_se@meta.data %>%
  tibble::rownames_to_column("barcodes") %>%
  dplyr::left_join(hist_df, by = c("barcodes", "gem_id")) %>%
  dplyr::mutate(
    histological_annotation = dplyr::if_else(is.na(histological_annotation),
      "unannotated", histological_annotation)) %>%
  tibble::column_to_rownames("barcodes")

table(merged_se@meta.data$histological_annotation)
## 
##    Mantle zones     unannotated Vascular spaces 
##             165           15971              88

Visualize the annotation

library(colorBlindness)

SpatialPlot(
  object = merged_se,
  group.by = "histological_annotation",
  images = c("esvq52_nluss5", "exvyh1_66caqq"),
  crop = FALSE, pt.size.factor = 1.25,
  cols = c("#009292", "lightgrey", "#ff6db6"))

Save object

"{anot}/{robj_dir}/integrated_spatial_annot.rds" %>%
  glue::glue() %>%
  here::here() %>%
  saveRDS(object = merged_se, file = .)

Session Info

sessionInfo()
## R version 4.0.1 (2020-06-06)
## Platform: x86_64-conda_cos6-linux-gnu (64-bit)
## Running under: Red Hat Enterprise Linux Server release 6.7 (Santiago)
## 
## Matrix products: default
## BLAS/LAPACK: /scratch/groups/hheyn/software/anaconda3/envs/spatial_r/lib/libopenblasp-r0.3.12.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=es_ES.UTF-8        LC_COLLATE=en_US.UTF-8     LC_MONETARY=es_ES.UTF-8    LC_MESSAGES=en_US.UTF-8    LC_PAPER=es_ES.UTF-8       LC_NAME=C                  LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=es_ES.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] colorBlindness_0.1.9 readr_2.1.1          stringr_1.4.0        dplyr_1.0.7          cowplot_1.1.1        ggpubr_0.4.0         ggplot2_3.3.5        SeuratObject_4.0.4   Seurat_4.0.5        
## 
## loaded via a namespace (and not attached):
##   [1] Rtsne_0.15            colorspace_2.0-2      ggsignif_0.6.3        deldir_1.0-6          ellipsis_0.3.2        ggridges_0.5.3        rprojroot_2.0.2       spatstat.data_2.1-0   farver_2.1.0          leiden_0.3.9          listenv_0.8.0         bit64_4.0.5           ggrepel_0.9.1         fansi_0.5.0           codetools_0.2-18      splines_4.0.1         knitr_1.36            polyclip_1.10-0       jsonlite_1.7.2        broom_0.7.10          ica_1.0-2             cluster_2.1.2         png_0.1-7             uwot_0.1.11           shiny_1.7.1           sctransform_0.3.2     spatstat.sparse_2.0-0 compiler_4.0.1        httr_1.4.2            backports_1.4.1       assertthat_0.2.1      Matrix_1.4-0          fastmap_1.1.0         lazyeval_0.2.2        cli_3.1.0             later_1.3.0           htmltools_0.5.2       tools_4.0.1           igraph_1.2.10         gtable_0.3.0          glue_1.5.1            RANN_2.6.1            reshape2_1.4.4        Rcpp_1.0.7            carData_3.0-4         scattermore_0.7       jquerylib_0.1.4       vctrs_0.3.8           nlme_3.1-153          lmtest_0.9-39         xfun_0.29             globals_0.14.0        mime_0.12             miniUI_0.1.1.1       
##  [55] lifecycle_1.0.1       irlba_2.3.5           rstatix_0.7.0         goftest_1.2-3         future_1.23.0         MASS_7.3-54           zoo_1.8-9             scales_1.1.1          vroom_1.5.7           spatstat.core_2.3-2   hms_1.1.1             promises_1.2.0.1      spatstat.utils_2.3-0  parallel_4.0.1        RColorBrewer_1.1-2    yaml_2.2.1            reticulate_1.22       pbapply_1.5-0         gridExtra_2.3         sass_0.4.0            rpart_4.1-15          stringi_1.7.6         highr_0.9             rlang_0.4.12          pkgconfig_2.0.3       matrixStats_0.61.0    evaluate_0.14         lattice_0.20-45       ROCR_1.0-11           purrr_0.3.4           tensor_1.5            labeling_0.4.2        patchwork_1.1.1       htmlwidgets_1.5.4     bit_4.0.4             tidyselect_1.1.1      here_1.0.1            parallelly_1.29.0     RcppAnnoy_0.0.19      plyr_1.8.6            magrittr_2.0.1        R6_2.5.1              generics_0.1.1        DBI_1.1.1             pillar_1.6.4          withr_2.4.3           mgcv_1.8-38           fitdistrplus_1.1-6    survival_3.2-13       abind_1.4-5           tibble_3.1.6          future.apply_1.8.1    car_3.0-12            crayon_1.4.2         
## [109] KernSmooth_2.23-20    utf8_1.2.2            spatstat.geom_2.3-1   plotly_4.10.0         tzdb_0.2.0            rmarkdown_2.11        grid_4.0.1            data.table_1.14.2     digest_0.6.29         xtable_1.8-4          tidyr_1.1.4           httpuv_1.6.4          gridGraphics_0.5-1    munsell_0.5.0         viridisLite_0.4.0     bslib_0.3.1